Section 15 Measure Classification

15.1 Introduction

Our measures are too loosely grouped, making it difficult for our optimization. We are aiming to have around 500 measures, currently we have more.

We will need our HRU ids spatially, as well as our Measures, erosion classes, and property register.

require(sf)

## Loading required package: sf

## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

require(dplyr)

## Loading required package: dplyr

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

require(mapview)

## Loading required package: mapview

require(DT)

## Loading required package: DT

Get our maps:

# For measure IDs and location:
lu_map <- read_sf("model_data/input/land/CS10_LU.shp")
# For HRU IDs:
hru_map <- read_sf("model_data/cs10_setup/optain-cs10/data/vector/hru.shp")
# For Erosion Risk Levels:
erosion_map <- read_sf("model_data/input/soil/erosionriskclasses_kraakstad.shp")
# For Farm ownership IDs:
farm_map <- read_sf("model_data/input/property/matrikkel.shp")
farm_map <- st_transform(farm_map, st_crs(lu_map))

st_crs(lu_map) == st_crs(hru_map)

## [1] TRUE

st_crs(erosion_map)==  st_crs(hru_map)

## [1] TRUE

st_crs(farm_map) ==  st_crs(hru_map)

## [1] TRUE

# filter data
lu_filter <- lu_map %>% select("buffer_6m_", "gully", "wetland", "dam", "type", "geometry")
hru_filter <-  hru_map  %>% select("name", "type")
erosion_filter <- erosion_map %>% select("A_HPKL", "geometry")
farm_filter <- farm_map %>% select("GARDSNUMME", "geometry")

farm_dissolved <- farm_filter %>%
  group_by(GARDSNUMME) %>%
  summarise()

erosion_dissolved <- erosion_filter %>%
  group_by(A_HPKL) %>%
  summarise()

mapview(farm_dissolved, zcol = "GARDSNUMME")

mapview(erosion_dissolved, zcol = "A_HPKL")

Ideally the above 4 maps would have their attributes spatially joined in R, but since I cannot figure out a nice way to do that, I did it quick and nicely in QGIS:

mea_map <- read_sf("model_data/input/measure/measure_comola_map.shp")
mapview(mea_map, zcol = "A_HPKL")

15.2 Determining Classification

We plan to do the following farm based classification:

70 farms in the catchment could implement gullies:

gully_farms <- mea_map %>% filter(gully > 0) %>% pull(GARDSNUMME) %>% unique() %>% length()
gully_farms

## [1] 70

68 farms can implement buffers:

buffer_farms <- mea_map %>% filter(buffer_6m_ > 0) %>% pull(GARDSNUMME) %>% unique() %>% length()
buffer_farms

## [1] 68

And we can implement 4 ponds.”

dams <- mea_map %>% pull(dam) %>% unique() %>% length()
dams

## [1] 4

We have 501 agricultural fields to apply no-till to.

non_agri <- c("frst", "past", "rngb", "urml", "utrn", "wetf")
mea_map %>% filter((type %in% non_agri) == FALSE) %>% pull(type) %>% unique() %>% length()

## [1] 501

If we group these by farm level, we have 82

mea_map %>% filter((type %in% non_agri) == FALSE) %>% pull(GARDSNUMME) %>% unique() %>% length()

## [1] 82

If we expand these to erosion risk, we have:

class1 = mea_map %>% filter((type %in% non_agri) == FALSE) %>% filter(A_HPKL == 1) %>% pull(GARDSNUMME) %>% unique() %>% length()

class2 = mea_map %>% filter((type %in% non_agri) == FALSE) %>% filter(A_HPKL == 2) %>% pull(GARDSNUMME) %>% unique() %>% length()

class3 = mea_map %>% filter((type %in% non_agri) == FALSE) %>% filter(A_HPKL == 3) %>% pull(GARDSNUMME) %>% unique() %>% length()

class4 = mea_map %>% filter((type %in% non_agri) == FALSE) %>% filter(A_HPKL == 4) %>% pull(GARDSNUMME) %>% unique() %>% length()

noclass = mea_map %>% filter((type %in% non_agri) == FALSE) %>% filter(is.na(A_HPKL)) %>% pull(GARDSNUMME) %>% unique() %>% length()

hrus_df <- data.frame (c1 = class1, c2 = class2, c3 = class3, c4 = class4, noclass=noclass)

DT::datatable(hrus_df)

  class1 + # low erosion risk lowtill farms
  class2 + # medium erosion risk lowtill farms
  class3 + # high erosion risk lowtill farms
  class4 + # very erosion risk lowtill farms
  noclass + # no erosion risk lowtill farms (not mapped)
  gully_farms + # gully farms
  buffer_farms + # buffer strip farms
  dams    # ponds

## [1] 364

15.3 Applying Classification

First we need to allocate incremental IDs to the measures, based

class1_alloc <- c(1:class1)
class2_alloc <- c((last(class1_alloc)+1):(last(class1_alloc)+class2))
class3_alloc <- c((last(class2_alloc)+1):(last(class2_alloc)+class3))
class_4_alloc <- c((last(class3_alloc)+1):(last(class3_alloc)+class4))
class_0_alloc <- c((last(class_4_alloc)+1):(last(class_4_alloc)+noclass))
gully_alloc <- c((last(class_0_alloc)+1):(last(class_0_alloc)+gully_farms))
buffer_alloc <- c((last(gully_alloc)+1):(last(gully_alloc)+buffer_farms))
dam_alloc <- c((last(buffer_alloc)+1):(last(buffer_alloc)+dams))

full_Series <- c(class1_alloc,class2_alloc,class3_alloc, class_4_alloc, class_0_alloc, gully_alloc, buffer_alloc, dam_alloc) 

full_Series %>% length() # same length?

## [1] 364

(full_Series-c(1:364) == 0) %>% all() # all incremental?

## [1] TRUE

now we need to translate the gardnummers to the allocated measure IDs..

15.3.1 Gullies (grassed water ways)

Measure range from 223 to 292

# Grabs the measure codes for each gully measure
gully_codes <- mea_map %>% filter(gully > 0) %>% pull(gully) 

# grabs the farm IDs for each gully measure
gully_farm_ids <- mea_map %>% filter(gully > 0) %>% pull(GARDSNUMME)

# Grabs the unique farms which have gully measures
gully_farm_uniq <- gully_farm_ids %>% unique() %>% sort()

# A look up table from gully farm IDs to old measure IDs
mea_to_farm <- data.frame(farm_id = gully_farm_ids, measure_id = gully_codes)

# A look up table from farm IDs which have gullies, to the gully-allocated IDs
farm_to_alloc <- data.frame(farm_id = gully_farm_uniq, gully_alloc = gully_alloc)

# A look up table combined from previous two
lookup <- left_join(farm_to_alloc,mea_to_farm, by = "farm_id", keep = FALSE, )

## Recoding IDs
# https://stackoverflow.com/questions/66231321/recode-values-based-on-look-up-table-with-dplyr-r

gully_df = mea_map %>% st_drop_geometry %>% filter(gully > 0) %>% select(id, name, GARDSNUMME, gully)

recode_vec <- lookup$gully_alloc
names(recode_vec) <- lookup$measure_id

# reclassifiying
gully_reclassed <- gully_df %>% mutate(gully_new = recode(gully, !!!recode_vec))

datatable(gully_reclassed)

rm(gully_df, gully_alloc, gully_codes, gully_farm_ids, gully_farm_uniq, gully_farms)

15.3.2 Buffer Strips

Same procedure, measure IDs 293 to 360

# Grabs the measure codes for each buffer_6m_  measure
buffer_codes <- mea_map %>% filter(buffer_6m_  > 0) %>% pull(buffer_6m_ ) 

# grabs the farm IDs for each buffer measure
buffer_farm_ids <- mea_map %>% filter(buffer_6m_  > 0) %>% pull(GARDSNUMME)

# Grabs the unique farms which have buffer_6m_  measures
buffer_farm_uniq <- buffer_farm_ids %>% unique() %>% sort()

# A look up table from buffer farm IDs to old measure IDs
mea_to_farm <- data.frame(farm_id = buffer_farm_ids, measure_id = buffer_codes)

# A look up table from farm IDs which have buffers, to the buffer-allocated IDs
farm_to_alloc <- data.frame(farm_id = buffer_farm_uniq, buffer_alloc = buffer_alloc)

# A look up table combined from previous two
lookup <- left_join(farm_to_alloc,mea_to_farm, by = "farm_id", keep = FALSE, )

## Recoding IDs
# https://stackoverflow.com/questions/66231321/recode-values-based-on-look-up-table-with-dplyr-r

buffer_df = mea_map %>% st_drop_geometry %>% filter(buffer_6m_ > 0) %>% select(id, name, GARDSNUMME, buffer_6m_)

recode_vec <- lookup$buffer_alloc
names(recode_vec) <- lookup$measure_id

# reclassifiying
buffer_reclassed <- buffer_df %>% mutate(buffer_new = recode(buffer_6m_, !!!recode_vec))

datatable(buffer_reclassed)

rm(buffer_df, buffer_alloc, buffer_codes, buffer_farm_ids, buffer_farm_uniq, buffer_farms)

15.3.3 NoTill (management measure)

# farm codes of farms with class 1 erosion agricultural fields.
mea_map %>% filter((type %in% non_agri) == FALSE) %>% filter(A_HPKL == 1) %>% select(name, GARDSNUMME, type, geometry)

## Simple feature collection with 636 features and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 604019.9 ymin: 6608002 xmax: 610154.4 ymax: 6622017
## Projected CRS: ETRS89 / UTM zone 32N
## # A tibble: 636 × 4
##    name    GARDSNUMME type                                              geometry
##    <chr>        <dbl> <chr>                                        <POLYGON [m]>
##  1 hru0009          1 a_001f_5 ((605305.1 6619679, 605307.2 6619677, 605309.8 6…
##  2 hru0010          3 a_002f_2 ((605503.8 6619824, 605508.2 6619824, 605512.2 6…
##  3 hru0011          3 a_002f_2 ((605390.2 6619491, 605336.8 6619456, 605326.8 6…
##  4 hru0014          3 a_002f_1 ((605400.2 6619337, 605414.5 6619303, 605425.2 6…
##  5 hru0053          3 a_002f_2 ((605404.8 6619515, 605404.5 6619515, 605392.7 6…
##  6 hru0055          3 a_002f_1 ((605440 6619243, 605439.9 6619242, 605439 66192…
##  7 hru0057          3 a_002f_2 ((605538.2 6619626, 605535 6619618, 605523.8 661…
##  8 hru0058          3 a_002f_2 ((605432.1 6619472, 605428.3 6619462, 605422.7 6…
##  9 hru0062          3 a_002f_1 ((605385.3 6619367, 605387.9 6619364, 605391.1 6…
## 10 hru0066          3 a_002f_2 ((605437.2 6619489, 605437.5 6619486, 605435.1 6…
## # ℹ 626 more rows

15.4 HRU Scenario Table

We now need to generate the scenario HRU table for COMOLA.

The land use shapefile must include columns for each measure (names according to the list above!) for each measure the id (as integer) of potential sites for implementation must be given.

Multiple hrus can have the same measure id, e.g. if an individual grassed waterway consists of multiple hrus if two individual measures of the same type (e.g. two crossing grassed waterways) occur on the same hru, provide both ids separated by ‘_’, e.g. ‘5_6’ (see also example data).

Define NA in case there is no potential site for a measure (please study the example lu shapefile)

# hru shapefile which was generated by SWATbuildR (.../data/vector/hru.shp)
hru_shp <- 'cs10_setup/optain-cs10/data/vector/hru.shp'

lu_shp <- ""
lu_shp_sf <- read_sf(lu_shp)

library(sf)
library(data.table)
library(stringr)
require(dplyr)


# list of all OPTAIN measures
measures_optain <- data.frame(
  'nswrm_code' = c(
    'buffer',
    'edgefilter',
    'hedge',
    'grassslope',
    'grassrchrg',
    'pond',
    'afforest',
    'floodres',
    'channres',
    'swale',
    'wetland',
    'cdrain',
    'terrace',
    'notill',
    'lowtill',
    'lowtillcc',
    'mulching',
    'subsoiling',
    'rotation',
    'intercrop',
    'covercrop',
    'earlysow',
    'droughtplt'
  ),
  'description' = c(
    'Riparian buffers',
    'Edge-of-field filter strips',
    'Hedges/Field division',
    'Grasland cover on erosive slopes',
    'Grasland cover in recharge area',
    'Retention/detention ponds',
    'Afforestation',
    'Floodplain restoration',
    'Channel restoration',
    'Swales',
    'Constructed wetlands',
    'Controlled drainage',
    'Terracing',
    'No-till agriculture',
    'Low-till agriculture',
    'Low-till agriculture combined with cover crops',
    'Mulching',
    'Subsoiling',
    'Crop rotation',
    'Intercropping',
    'Green cover/ catch crops',
    'Early sowing',
    'Drought-resistant plants'
  )
)

lu_shp_sf = rename(lu_shp_sf, buffer = buffer_6m_)
lu_shp_sf$buffer[which(lu_shp_sf$buffer == 0)] = NA

lu_shp_sf = rename(lu_shp_sf, grassslope = gully)
lu_shp_sf$grassslope[which(lu_shp_sf$grassslope == 0)] = NA

lu_shp_sf = rename(lu_shp_sf, pond = dam)
lu_shp_sf$grassslope[which(lu_shp_sf$grassslope == 0)] = NA


lu_shp_sf$wetland[which(lu_shp_sf$wetland == 0)] = NA

# generate scenario hru table (just run the code)--------------------------

is.integer0 <- function(x)
{
  is.integer(x) && length(x) == 0L
}

lu <- lu_shp_sf
hru_swat <- read_sf(hru_shp)

# assign measures to SWAT+ hru names
lu_swat <- st_drop_geometry(st_join(hru_swat, lu, largest=T))
hru <- data.frame('hru'=lu_swat$name)

# generate hru-scenario table
meas.col <- which(colnames(lu_swat) %in% measures_optain$nswrm_code)
hru_scn <- data.frame('hru'=NA, 'scenario_id'=NA)

for(i in 1:length(meas.col)){
  idx <- which(!is.na(lu_swat[,meas.col[i]]))
  meas_ <- colnames(lu_swat)[meas.col[i]]
  hru_ <- hru
  hru_$scenario_id <- NA
  hru_$scenario_id[idx] <- paste0(meas_,'_', as.matrix(lu_swat[idx,meas.col[i]]))
  # different grassed waterways on the same hru?
  multiple <- which(str_count(hru_$scenario_id, "_") > 1)
  if(!is.integer0(multiple)){
    scn_multiple <- str_split_fixed(hru_$scenario_id[multiple],"_",3)
    hru_$scenario_id[multiple] <- paste0(meas_,'_',scn_multiple[,2])
    hru_add <- data.frame('hru' = hru_$hru[multiple],
                          'scenario_id' = paste0(meas_,'_',scn_multiple[,3]))
    hru_ <- rbind.data.frame(hru_, hru_add)
  }
  hru_ <- hru_[-which(is.na(hru_$scenario_id)),]
  hru_scn <- rbind.data.frame(hru_scn, hru_)
}

hru_scn <- hru_scn[-1,]
rownames(hru_scn) <- c(1:dim(hru_scn)[1])
ind.list <- unique(hru_scn$scenario_id)
measures <- setNames(replicate(length(ind.list),data.frame()),ind.list)

for(m in 1:length(ind.list)){
  measures[[m]] <- hru_scn$hru[which(hru_scn$scenario_id==ind.list[m])]
}

measures_df <- data.frame(lapply(measures, function(x) {
  x <- unlist(x)
  length(x) <- max(lengths(measures))
  return(x)
}))

measures_df_t <- data.table::transpose(measures_df)
measures_df_t$type <- colnames(measures_df)
measures_df_t <- measures_df_t[,c(length(measures_df_t),1:(length(measures_df_t)-1))]

# write hru-scenario table
write.csv(measures_df_t, 'scenario_hru_CS10.csv',quote=F, row.names=F)